DOUBLE PRECISION function y_Euler(x0, y0, b, h, fungsi)
	IMPLICIT NONE
	double precision x0, y0, b, h, fungsi
	double precision x, y
	integer*2 r, n	
	n = (b-x0)/h
	y = y0
	x = x0
	r = 1
	while (r<=n) do
		y = y + h*fungsi(x,y)
		x = x + h
		r = r + 1
	end do
	y_Euler = y
end


DOUBLE PRECISION function y_Heun(x0, y0, b, h, fungsi)
	IMPLICIT NONE
	double precision x0, y0, b, h, fungsi
	double precision x, y, y_s
	integer*2 r, n
	n = (b-x0)/h
	y = y0
	x = x0
	r = 1
	while (r<=n) do
		y_s = y
		y = y + h*fungsi(x,y)
		y = y_s + (h/2)*(fungsi(x,y_s)+fungsi(x+h,y))
		x = x + h
		r = r + 1
	end do
	y_Heun = y
end

subroutine y_Euler_param3(x0,y0,z0,b,h,fungsi_y,fungsi_z,y_hasil,z_hasil)
	IMPLICIT NONE
	double precision x0, y0, z0, b, h, fungsi_y, fungsi_z
	double precision x, y, z, y_hasil, z_hasil, temp
	integer*2 r, n	
	n = (b-x0)/h
	y = y0
	x = x0
	z = z0
	r = 1
	while (r<=n) do
		temp = y
		y = y + h*fungsi_y(x,y,z)
		z = z + h*fungsi_z(x,temp,z)
		x = x + h
		r = r + 1
	end do
	y_hasil = y
	z_hasil = z
end


subroutine y_Heun_param3(x0,y0,z0,b,h,fungsi_y,fungsi_z,y_hasil,z_hasil)
	IMPLICIT NONE
	double precision x0, y0, z0, b, h, fungsi_y,fungsi_z
	double precision x, y, z, y_s, z_s, y_hasil, z_hasil
	integer*2 r, n
	n = (b-x0)/h
	y = y0
	x = x0
	z = z0
	r = 1
	while (r<=n) do
		y_s = y
		y = y + h*fungsi_y(x,y,z)
		y = y_s + (h/2)*(fungsi_y(x,y_s,z)+fungsi_y(x+h,y,z))
		z_s = z
		z = z + h*fungsi_z(x,y_s,z)
		z = z_s + (h/2)*(fungsi_z(x,y_s,z_s)+fungsi_z(x+h,y_s,z))
		x = x + h
		r = r + 1
	end do
	y_hasil = y
	z_hasil = z
end


subroutine runge_kutta_param3(x0,y0,z0,x_tanya,jarak,fungsi_y,fungsi_z,y_hasil,z_hasil)
    IMPLICIT NONE
	double precision x0, y0, z0, x_tanya, jarak,fungsi_y,fungsi_z
	integer*2 r,n
	double precision k1,k2,k3,y_hasil,z_hasil,temp,x,y,z
	n = (x_tanya-x0)/jarak
	r=1
	x = x0
	y = y0
	z = z0
	while (r<=n) do
		temp = y
		k1 = jarak*fungsi_y(x,y,z)
		k2 = jarak*fungsi_y(x+jarak/2,y+k1/2,z+k1/2)
		k3 = jarak*fungsi_y(x+jarak,y-k1+2*k2,z-k1+2*k2)
		y = y + (k1+4*k2+k3)/6
		k1 = jarak*fungsi_z(x,temp,z)
		k2 = jarak*fungsi_z(x+jarak/2,temp+k1/2,z+k1/2)
		k3 = jarak*fungsi_z(x+jarak,temp-k1+2*k2,z-k1+2*k2)
		z = z + (k1+4*k2+k3)/6
		x = x + jarak
		r=r+1
	end do
	y_hasil = y
	z_hasil = z
end

subroutine runge_kutta_param4(x0,y0,z0,x_tanya,jarak,fungsi_y,fungsi_z,y_hasil,z_hasil)
    IMPLICIT NONE
	double precision x0, y0, z0, x_tanya, jarak,fungsi_y,fungsi_z
	integer*2 r,n
	double precision k1,k2,k3,k4,y_hasil,z_hasil,temp,x,y,z
	n = (x_tanya-x0)/jarak
	r=1
	x = x0
	y = y0
	z = z0
	while (r<=n) do
		temp = y
		k1 = jarak*fungsi_y(x,y,z)
		k2 = jarak*fungsi_y(x+jarak/2,y+k1/2,z+k1/2)
		k3 = jarak*fungsi_y(x+jarak/2,y+k2/2,z+k2/2)
		k4 = jarak*fungsi_y(x+jarak,y+k3,z+k3)
		y = y + (k1+2*k2+2*k3+k4)/6
		k1 = jarak*fungsi_z(x,temp,z)
		k2 = jarak*fungsi_z(x+jarak/2,temp+k1/2,z+k1/2)
		k3 = jarak*fungsi_z(x+jarak/2,temp+k2/2,z+k2/2)
		k4 = jarak*fungsi_z(x+jarak,temp+k3,z+k3)
		z = z + (k1+2*k2+2*k3+k4)/6
		x = x + jarak
		r=r+1
	end do
	y_hasil = y
	z_hasil = z
end

subroutine adamsmuolton_param3(x0,y0,z0,x_tanya,h,fungsi_y,fungsi_z,y_hasil,z_hasil)
    IMPLICIT NONE
	double precision x0, y0, z0, x_tanya,h,fungsi_y,fungsi_z
	integer*2 r,n
	double precision y_hasil,z_hasil,temp,temp2,x,y,z
	double precision y1,y2,y3,z1,z2,z3,p1,p2
	external runge_kutta_param3
	n = (x_tanya-x0)/h
	if (n<4) then
		print*,"N = ",n," kurang dari 4."
		y = 9999
		z = 9999
	else
		r=1
		x = x0
		y = y0
		z = z0
		call runge_kutta_param3(x,y,z,x+h,h,fungsi_y,fungsi_z,y1,z1)
		call runge_kutta_param3(x,y,z,x+2*h,h,fungsi_y,fungsi_z,y2,z2)
		call runge_kutta_param3(x,y,z,x+3*h,h,fungsi_y,fungsi_z,y3,z3)
		x = x+3*h
		r = 4
		while (r<=n) do
			p1=y3+h/24*(-9*fungsi_z(x-3*h,y,z)+37*fungsi_y(x-2*h,y1,z1)-59*fungsi_y(x-h,y2,z2)+55*fungsi_y(x,y3,z3))
			p2=z3+h/24*(-9*fungsi_z(x-3*h,y,z)+37*fungsi_z(x-2*h,y1,z1)-59*fungsi_z(x-h,y2,z2)+55*fungsi_z(x,y3,z3))
			temp = p1
			temp2 = p2
			p1=y3+h/24*(fungsi_y(x-2*h,y1,z1)-5*fungsi_y(x-h,y2,z2)+19*fungsi_y(x,y3,z3)+9*fungsi_y(x+h,temp,temp2))
			p2=z3+h/24*(fungsi_z(x-2*h,y1,z1)-5*fungsi_z(x-h,y2,z2)+19*fungsi_z(x,y3,z3)+9*fungsi_z(x+h,temp,temp2))
			y = y1
			y1 = y2
			y2 = y3
			y3 = p1
			z = z1
			z1 = z2
			z2 = z3
			z3 = p2
			x = x + h
			r = r + 1
		end do
		y = p1
		z = p2
	end if
	y_hasil = y
	z_hasil = z
end


SUBROUTINE RUNGE_KUTTA_ORDE_3(f,x0, y0, b, h,output)
	IMPLICIT NONE
	DOUBLE PRECISION f,x0, y0, b, h, output, x, y, k1, k2, k3
	INTEGER r, n
	
	n = (b - x0)/h
	y = y0
	x = x0
	do r=1,n
		k1 = h*f(x,y)
		k2 = h*f(x+h/2,y+k1/2)
		k3 = h*f(x+h, y-k1+2*k2)
		y = y + (k1+4*k2+k3)/6
		x = x+h
	end do
	output = y;
END

SUBROUTINE RUNGE_KUTTA_ORDE_4(f,x0, y0, b, h,output)
	IMPLICIT NONE
	DOUBLE PRECISION f,x0, y0, b, h, output, x, y, k1, k2, k3, k4
	INTEGER r, n
	
	n = (b - x0)/h
	y = y0
	x = x0
	do r=1,n
		k1 = h*f(x,y)
		k2 = h*f(x+h/2,y+k1/2)
		k3 = h*f(x+h/2, y+k2/2)
		k4 = h*f(x+h, y+k3)
		y = y + (k1+2*k2+2*k3+ k4)/6
		x = x+h
	end do
	output = y;
END

SUBROUTINE ADAM_BASHFORTH_MOULTON(f,x0,y0,b,h,output)
	IMPLICIT NONE
	DOUBLE PRECISION x0,y0,b,h,output,x,y,y1,y2,y3,f
	INTEGER r,n
	EXTERNAL RUNGE_KUTTA_ORDE_3
	
	n = (b-x0)/h
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+h, h,y1)
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+2*h, h,y2)
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+3*h, h,y3)
	x = x0 + 3*h
	do r=4,n
		y = y3 + h/24*(-9*f(x-3,y0) + 37*f(x-2*h,y1) - 59*f(x-h,y2) + 55*f(x,y3))
		y = y3 + h/24*(f(x-2*h,y1) - 5*f(x-h,y2) + 19*f(x,y3) + 9*f(x+h,y))
		y0 = y1
		y1 = y2
		y2 = y3
		y3 = y
		x = x+h
	end do
	output = y
END